Quantum Fluctuations in the Superfluid State of the BCS-BEC Crossover 
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We determine the effects of quantum fluctuations about the T = mean field solution of the 
BCS-BEC crossover in a dilute Fermi gas using the functional integral method. These fluctuations 
are described in terms of the zero point motion of collective modes and the virtual scattering of 
gapped quasiparticles. We calculate their effects on various measurable properties, including chem- 
ical potential, ground state energy, the gap, the speed of sound and the Landau critical velocity. 
At unitarity, we find excellent agreement with quantum Monte Carlo and experimental results. In 
the BCS limit, we show analytically that we obtain Fermi liquid interaction corrections to ther- 
modynamics including the Hartree shift. In the BEC limit, we show that the theory leads to an 
approximate description of the reduction of the scattering length for bosonic molecules and also 
obtain quantum depletion of the Lee- Yang form. At the end of the paper, we describe a method 
to include feedback of quantum fluctuations into the gap equation, and discuss the problems of 
self-consistent calculations in satisfying Goldstone's theorem and obtaining ultraviolet finite results 
at unitarity. 
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I. INTRODUCTION 

The BCS-BEC crossover [l[3,i,ll is a problem of long 
standing interest in many-body physics with implications 
for a variety of fields including condensed matter, high- 
energy, nuclear and atomic and molecular physics. Re- 
cent experimental progress in cooling atomic Fermi gases 
to ultralow temperatures and tuning the interactions be- 
tween atoms using the Feshbach resonance technique has 
led to an explosion of interest in the BCS-BEC crossover 

SiSiiEEi. 

The theoretical problem is to determine the properties 
of a system with two species of fermions, spin-up and 
down, with equal masses and densities, interacting via a 
short range attractive potential described by a scatter- 
ing length as- The two extremes of the crossover are 
well understood theoretically. Weak attractive interac- 
tions characterized by a small, negative scattering length 
as lead to collective Cooper pairing of atoms and BCS 
superfluidity. In the opposite limit of large attraction, 
characterized by a small, positive scattering length as 
one obtains bosonic molecules which exhibit BEC. The 
intermediate regime, around the unitary point at which 
I as I — > oo, where one has a strongly interacting Fermi 
gas, is the most interesting and least well understood 
theoretically. 

The original mean field (MF) theory of Leggett Q and 
Eagles [2] does a decent job of describing the entire T ~ 
crossover at a qualitative level [l^ . It has only one addi- 
tional ingredient to the standard BCS theory: the chem- 
ical potential must be determined self-consistently along 
with the pairing gap. This is sufficient to give qual- 
itatively reasonable results 14| which evolve smoothly 
through unitarity all the way up to the molecular BEC. 

Recent theoretical and experimental developments 
have led to a realization of the quantitative shortcomings 
of mean field (MF) theory at T = especially at unitar- 
ity. The ground state energy density at unitarity is of 



the form £qJN = (1 -|- /3)(3e//5), which is a "universal" 
number (rl.flit times the free Fermi gas energy since there 
is no scale other than ef as |a g| — > oo. Quantum Monte 
Carlo (QMC) calculations ^,[l3l obtain {l+p)= 0.44, 
while experiments [8, 9, lOj find (1 -I- P) in the range 
0.32 to 0.44. In contrast the ground state energy density 
within MF theory yields {I + P) = 0.59, which is 
about 34% larger than the QMC result. Furthermore, in 
the BEC hmit [kpas 0+), although the MET correctly 
predicts a repulsive interaction between the constituent 
bosons, it misses 0{kpas)^^^ corrections to thermody- 
namic quantities , which are present in a weakly repulsive 
Bose gas. 

Our main motivations were to understand this quan- 
titative discrepancy, on which there has been recent 
progress by several approaches 0, 0, El, S lii (dis- 
cussed below) , and also to get a physical picture of the 
quantum fluctuations missing in MF theory that are re- 
sponsible for such a large energy difference. We show here 
that the many-body ground state in the crossover must 
include, in addition to BCS pairing, the effects of the zero 
point motion of the collective excitations - the oscillation 
of the phase (and amplitude) of the order parameter - 
and the effects of virtual scattering of quasiparticle exci- 
tations. 

Our central result, from which essentially all our other 
results follow, is that the thermodynamic potential at 
T = is given by 
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Here the first line represents the "fermionic contribu- 
tion" to the ground state energy of the superfluid. It 
has the same structure as the mean field result and may 
be thought of as coming from filling up the negative en- 
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Mean Field 
Theory [141 


0.59 


0.44 


0.69 


Quantum 
Monte Carlo [17, M 


0.44 


0.38 


0.6 ±0.1 


Experiments 
[8, 9, 


0.32 - 0.51 


0.38 




Mean Field + 
Fluctuations 


0.40 


0.37 


0.47 



TABLE I: Comparison of ground state energy, speed of sound 
and gap at unitarity obtained by different methods. The last 
row gives the results obtained in sections V and Vlll. The 
results of a self-consistent calculation are described in section 
IX. 



ergy states in a Bogoliubov deGennes framework. The 
second line represents the "bosonic contribution" which 
arises from Gaussian fluctuations about the saddle point. 
It consists of the zero-point energy of the collective mode 
with dispersion wo(q) and of an integral which describes 
the contribution from the virtual scattering of quasiparti- 
clcs (with a phase shift S{q,uj)), whose two-particle con- 
tinuum begins at i?c(q). The last term TZ regularizes 
the ultraviolet divergence in the bosonic contribution and 
will be described in detail later; see eq. (|27p . 

The problem of determining the ground state energy 
density of a strongly interacting system is analogous to 
that of determining the cosmological constant in quan- 
tum field theories. The latter problem is notorious for 
being dominated by physics at the scale of the ultraviolet 
cutoff. Here we will show that our results are indepen- 
dent of the momentum cutoff, which is the inverse of the 
range roff of the attractive potential between fermions. 

Although there has been considerable attention de- 
voted to "universality" at unitarity {\as\ ~ oo), we em- 
phasize the simple but often overlooked point that any 
observable quantity at T = is a universal function of the 
single parameter 1/kpas- By universal we mean that the 
results are independent of microscopic details below the 
ultraviolet cutoff length scale of r^s, provided kpreS <Si 1. 
Thus it does not matter whether one looks at an exper- 
iment with ^Li or ""^K, the result can only depend on 
system parameters through the combination l/kpas- In 
the absence of a small parameter in the crossover prob- 
lem, we judge the validity of our approximations not only 
by their success at unitarity in comparison with quantum 
Monte Carlo and experiments, but also by their ability 
to reproduce known results in the BCS and BEG limits. 

We conclude this Section with a summary of our main 
results and an outline of the rest of the paper. In Sec- 



tions II, III and IV we describe the functional integral 
formalism used in the paper and results for the chemical 
potential, gap, ground state energy, speed of sound and 
Landau critical velocity function of l/{kpas) are 
presented in Section V. Our main results include: 

(a) In the extreme BCS limit, the fluctuation correc- 
tions are dominated by the virtual scattering of ferniionic 
quasiparticles. We show in Section VI that we recover the 
exact Fermi-liquid corrections to the thermodynamics of 
a dilute Fermi gas, which are the Hartree shift of order 
kpo-s and the Galitskii and Huang-Lee- Yang corrections 



2j,|25|, |26| of order {kpas)"^, albeit with a negative scat- 
tering length as- We note that these are obtained from 
Gaussian fluctuations about the broken symmetry state 
and not by including them in an ad hoc way. 

(b) In the extreme BEC limit, the zero-point motion 
of the Bogoliubov sound mode dominates the thermody- 
namics. From the leading order corrections we estimate 
in Section VII the effective scattering length between the 
molecular bosons to be ~ 0.55as, an approximate re- 
sult which turns out to be close to the exact result 

for the four-body problem of O.Gag. At the next order 
we recover the Lee- Yang form for the quantum depletion 
[25I . [2^ of the molecular Bose gas with a coefficient that 
is only 6% less than the correct asymptotic expression. 

(c) At unitarity, both the zero-point motion of the 
collective modes and virtual scattering of quasiparticles 
are important. Our numerical results (see Table I) for the 
ground state energy, the gap and the speed of sound are 
in good agreement with experimental data and Quantum 
Monte Carlo results; See Section VIII for details. 

(d) The critical velocity Vc across the crossover is max- 
imum near unitarity, as previously predicted H, 2§| and 
as has been observed in experiments [30| . We estimate an 
upper bound on Vc using the Landau criterion, and find 
that quantum fluctuations considerably lower it with re- 
spect to mean field values; see Fig. 4. 

(e) The results described above are obtained within a 
scheme in which the Gaussian fluctuations do not feed 
back into the saddle point equation for the functional 
integral, and they only contribute to the thermodynamic 
potential ([T]). This is a natural approximation within the 
functional integral framework, and we show in Appendix 
G that it leads to exact answers in the simpler problem 
of the dilute repulsive Bose gas. 

(f) To go beyond this approximation, we next include 
in Section IX the self-consistent feedback of the Gaussian 
fluctuations into the gap equation. We find that this ap- 
proach leads to several problems, some of which we can 
resolve. For instance, we show how the apparent viola- 
tion of Goldstone's theorem in the self-consistent scheme 
can be resolved by going to an amplitude-phase represen- 
tation of the fluctuations. However, we point out that 
there are other problems which are not under control, 
such as the instability of the system in the extreme BEC 
limit. Our detailed analysis of the theory with a gap 
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equation modified by Gaussian corrections shows that 
imposing self-consistency does not necessarily lead to an 
improved approximation scheme. 

In Section X we compare our approach and results 
with several other methods which have been used to at- 
tack the same problem. Our ap proach has similarities 
with the 1 /N expansion [2l|, [24I but there are also dif- 
ferences which are discussed in Sec. X. The equations 
solved in Section V are the same as that obtained from 
the diagrammatic approach of Hu, Liu and Drummond 
[l^ . however our derivation is different and shows why 
it is natural not to renormalize the saddle-point condi- 
tion with fluctuation corrections. Further our approach 
also allows us to see what the impact of going beyond 
this approximation is, as indicated in (f) above. This 
gives insights into problems faced in other self-consistent 
calculations [23| . Our main conclusions are summarized 
at the end in Section XI. Technical details of the calcu- 
lations presented in the text are given in a series of six 
Appendices. In a seventh Appendix we use illustrate the 
methods used in the text for paired Fermi superfluids for 
the simpler case of a Bose superfluid. 



II. FUNCTIONAL INTEGRAL FORMALISM 

We consider a system of fermions of two species, which 
we call "spin" a ~^,[, each of mass m, described by the 
Hamiltonian density 



2™ 



(2) 



The first term has an implicit sum on the repeated in- 
dex cr, and the chemical potential fi is tuned to fix the 
average particle density n — fc|./(37r^) in a unit volume. 
Throughout the paper, we set h = ks = ^■ 

We consider the experimentally relevant case of a 
"broad" Feshbach resonance which can be adequately de- 
scribed within a single-channel formulation of a dilute 
gas with kpr^ff <C 1 where r^ff is the range of the poten- 
tial [sij . The two-body interaction in ([2|) is described by 
a "bare" coupling constant g and a momentum cut-off 
A, not explicitly shown above, which is of the order of 
l/rcff. The effective interaction at low-energies is com- 
pletely described by the s-wave scattering length as for 
the two-body problem in vacuum. To obtain a given 
renormalized as, the bare coupling 5(A) must be tuned 
using the relation 



Anaf 



-1 
5(A) 



E- 



|k|<A 



(3) 



where ek = |kp/2m. We will write most of our equa- 
tions in terms of the bare g, and only at the end we will 
use eq. to take the A ^ 00 limit and express the final 



results in terms of as- (For a detailed discussion of jus- 
tifying this regularization procedure, we refer the reader 
to Sec. IV of ref. [H). 

The partition function Z in the grand canonical ensem- 
ble at temperature T, chemical potential /x and in a unit 
volume, can be written as the imaginary time functional 
integral 3, 14 , 1^ over the Grassman fields ip and ip 



Z = j D4'^Di', exp {-3.4,) 
with the action 

5"^= dx {i;^{x)dripa{x) + H[tlj,i'] 



(4) 



(5) 



We use X to denote x = (x, r) where x is the spatial 
coordinate and r is imaginary time in the interval < 
T < (3, where /3 = l/T. J dx ^ dr J d'^x denotes an 
integral over all space and over imaginary time. Even 
though we are eventually interested in T = 0, we find it 
convenient to use the finite T Matsubara formalism and 
take /3 — > 00 at the end. 

We next use a Hubbard Stratanovich transformation 
with an auxiliary field A(a;) which couples to ip-^{x)'ip^{x) 
to obtain 

Z 



Dij^D^.DADA* exp (-^^.a) • 



(6) 



Using the spinor ip''{x) = (■0|(x), ■0|(a;)) and its hermi- 
tian conjugate ipi^) the action can be written as 



S. 



tp.A 



dx 



\A{x) 
9 



dxdx'i;\x)G-\x,x')ilj{x') 



(7) 

where the inverse Nambu-Gorkov Green's function G ^ 
is given by 



A*ix) 



A(x) 
-Or - V2/2m - fi 



X 6{x — x'). 

The functional integral is now quadratic in the fermion 
fields and these can be integrated out to obtain 



Z = J DAD A* exp (-^a) 



(9) 



with the action 



5a 



dx 



\A(x)l 
9 



TrlnG-i[A(a:)] (10) 



where the trace is over two-dimensional Nambu space. 
Eq. ([5]) is a formally exact expression for Z — exp(— /3r2), 
where 51 is the thermodynamic potential. 

III. MEAN FIELD THEORY 

We briefly discuss the mean fleld theory P, [l3| of the 
BCS-BEC crossover to introduce notation used through- 
out the paper. Technical details highlighting aspects 
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(such as convergence factors) which will be useful later 
are given in Appendix A. 

We begin by finding a spatially uniform, static saddle 
point Ao to the functional integral defined by eqs. (|9ll0p 
This is determined by the gap equation 



SSa/SAq = 0, 



(11) 



where 

5a [Ao] = 

with 

G-i[Ao] = 



^ TrlnGoi(fc) = 5o (12) 



Ao 



Ao 



^Go'ik)- (13) 



Here ikn ~ (2n+l)7r//3 are fermionic Matsubara frequen- 
cies, and = Ck — M with ek = |kp/2TO. 

After some straightforward algebra (see Appendix A) 
the T = gap equation (fTT|) can be finally written as 



= E 



1 



1 



2ek 2i?k 



(14) 



where = \/^^ + Aq and we have used eq. ^ to elim- 
inate g in favor of as- 

To determine both Aq and the chemical potential /i 
we need to use n = —dVl/dfi, in addition to (|14p . At the 
level of the mean field (MF) approximation, the thermo- 
dynamic potential is given by its saddle point estimate 
fio = Sq/ (3 which leads to the T = MF number equa- 
tion 



E 



1 



in 
^k 



(15) 



Equations pi]) and ([TS]) are the Leggett mean field equa- 
tions [H for the T = BCS-BEC crossover which can 
be solved to obtain the mean field values Ao and /i as a 
fimction of {kpas)' 



Finally, we can obtain an explicit result for the MF 
thermodynamic potential at T = in terms of Aq and 
/i. We evaluate the Matsubara sum in eq. (fT2|) . take the 
T = limit, and use eq. ([3]) to obtain 



r!o = 



IV. GAUSSIAN FLUCTUATIONS 

To go beyond the MF approximation and include the 
effects of fluctuations, we write 



A(a;) = Ao + ??(a;) 



(17) 



where the complex bosonic field field ri{x) describes 
space-time dependent fluctuations about the real, (x, r)- 
independent saddle point Aq . We Fourier transform from 
X = (x, r) — *■ g = (q, iqi) where iqi = i2-Kl/ (5 is the Mat- 
subara frequency for the bosonic 77 fields. We then write 



and 



as G ^ = Gq ^ + K, where Gg ^ is defined in H 



¥.{k,k + q) 











(18) 



We next expand the action Sa to order 77^. The first 
order term vanishes by the saddle point condition pT|) 
and we obtain 



Sa — Sq + Sg + . . . 

where the mcan-fleld Sq was deflned in eq. ([1 
Gaussian piece has the form 



ri*{~q) 



(19) 
The 

(20) 



The inverse fluctuation propagator M is given by [3, [s^ 
Mn(g) = M22(-9) = ^ + E G°2(fc)G°i(fc + q) 



E + E' 



(21) 



and 



Mi2(g) = M2i(<z) = J2 G?2(fc)G?2(fc + 

1 1 



E^ 



iqi+E + E' iqi- E- E' 



Here we use standard BCS notation 



vi = 1 - ui = I a 



(22) 



(23) 



together with the abbreviations u = u\^,v = v\^, E = 
Ek and u' = Mk+q, = fk+q, E' = -Ek+q- The flrst 
line in (PT|) and is valid at all temperatures, and 
the Matsubara sums lead to expressions involving Fermi 
hmctions (1 - / - /') and (/ - /'); see ref.Q- In the 
second line of ([2T|l and ((22|) we only give results valid in 
the r = hmit where both f{E) = f{E') = 0. The 
factor of 1 /g in (PT|) has to be regularized as usual using 
eq. ©. 

Integrating out the Gaussian fluctuations in 



Z ~ exp {-So) j DriD'n^ exp {-Sg) 



(24) 



we obtain an improved estimate of the thermodynamic 
potential 



"o + ^Eli^DetM(<z). 



(25) 
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FIG. 1: Spectrum of excitations which contribute to the Gaus- 
sian correction to the thermodynamic potential (|27p . The re- 
sults shown correspond to unitarity as ~ oo with fi = 0.4 ep 
and Ao — 0.465 ef- The full line shows the collective mode 
dispersion ajo(g) (pole of 1/DetM) and the shaded region is 
the two-particle continuum (branch cut). The dashed lines 
uj22{q) and Wii(g) are the zeros of , M22, and Mn respec- 
tively. 



where fig was defined in (|16p . We note that this result 
is true even for a non-Herniitian matrix M, provided its 
Hermitian part is positive definite 3J|. In our case, this 
condition corresponds to M11-I-M22 — 2Mi2 > 0, which is 
true whenever p4|) is satisfied. Physically this is related 
to an increase in energy under a distortion of the phase, 
as can be seen from the analysis of Section IX. 

There is however a problem with this expression (j25p . 
since it is actually ill-defined: the Matsubara sum is di- 
vergent and we need appropriate convergence factors to 
make it meaningful as discussed in detail in Appendix B. 
We only write the final result here: 



1 

2j3 



Mn(g) 
M22(g) 



DetM(q) 



(26) 



In order to gain physical insight into what eq. (|26p 
means, we will analytically continue from Matsubara fre- 
quencies to real frequencies: iqi ^ co + iO^ . Using stan- 
dard manipulations (see Appendix B) the Gaussian part 
of the thermodynamic potential at T = can be writ- 
ten as fig = -l/2J2qjl^duj/TT[S{q,u;) - (522(q,w) + 
(5ii(q, w)]. Here 6 is the phase of DetM defined by 
(5(q, w) = InilnDctM(q, w + iO+) and S22 and Su are 
the phases of M22 and Mn respectively. The integral 
runs only over oj < because at T = the Bosc factor 

The analytical structure of DetM(q, z) is as follows: It 
has zeros on the real axis dX z = ±a;o(q), which corre- 
spond to poles of the fiuctuation propagator, and describe 
the spectrum of collective excitations. These excitations 
are oscillations of the phase of the order parameter as 



5^0, and is the Goldstone mode arising from the bro- 
ken symmetry in the superfiuid state. We will show that 
^0(9) = as q ^ Q characteristic of a sound mode. In 
addition, at higher energies there are branch cuts along 
the real axis at each q, with branch points at ±i5c(q) 
with i?c(q) = min(iJk + -Ek+q)- These branch cuts rep- 
resent the two-particle continuum of states for scattering 
of gapped quasiparticles; see the lower panel of Fig. 11 
in Appendix D. 

On the negative u) axis, M22(q, w) and Mii(q, w) have 
zeros at — <x'22(q) and — a;ii(q) respectively, and each has 
its own scattering continuum. Although the physical 
meaning of these quantities is less clear, the role that they 
play in cutting off the ultraviolet divergences in eq. (|26p 
will be clarified in detail below. 

To illustrate these ideas, we show in Fig. [1] the col- 
lective mode spectra and the two-particle continuum at 
the unitary point at which \as\ diverges. Note that the 
collective mode frequency is initially linear in q, 

as expected, while the frequencies ti;22 and loh have non- 
zero values in the limit q Q. All of these frequencies 
eventually hit the two-particle continuum. Although we 
keep the integral in the Gaussian part fig over cj < 0, as 
it appears in the algebra, we find it simpler to plot all 
spectra as for positive excitation energies in Fig. [T] and 
subsequent figures. 

We note in passing that even though the unitary Fermi 
gas is a very strongly interacting system, nevertheless its 
collective mode spectrum does not show a roton-like min- 
imum observed in superfiuid Helium 4. We can under- 
stand this within a Feynman approach where the roton 
minimum arises from a peak in the static structure factor 
characteristic of a liquid, while here we are dealing with 
a gas, even if it is a very strongly interacting gas. 

Next we explicitly separate out the collective mode and 
continuum contributions and write the thermodynamic 
potential fi(r = 0) = £ - /iiV as 

fi = fio + i^[cJo(q)-t^22(q)+c^ii(q)-£^c(q)] (27) 
q 

~ 2^Ey dw[5(q,w)-(522(q,w)+(5ii(q,w)]. 



This is the full expression for the result ((T|) in the Intro- 
duction. 

The various contributions to the thermodynamic po- 
tential (|27p are now much more transparent compared 
with the Matsubara axis expression (|26p . fio is the mean 
field contribution to the ground state energy £. It 
may be interpreted as arising from filling up the nega- 
tive energy (— -Ek) fermionic states of the BCS Hamil- 
tonian, as is made clear in Bogoliubov-deGennes theory. 
The Gaussian contribution to £ has three parts to it. 
The first part a;o(q)/2 comes from the zero-point mo- 
tion of the collective mode. The second part, related 
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FIG. 2: Gap Aq and chemical potential ^ as a function of 
— {kpas)'^ across the BCS-BEC crossover. The dashed line 
is the mean field solution while the results of the calculation 
which includes Gaussian fluctuations are shown as solid lines. 

to the (5(q, w) terms, arises from virtual scattering of 
the fermionic quasiparticles whose two-particle contin- 
uum begins at the energy Ec{q). The third set of contri- 
butions, related to the W22,wii,(522 and terms, come 
from the convergence factors of eq. ((26|) and are essential 
to get a finite answer for fl. 

To get a better feel for these various contributions it is 
useful to look at limiting cases. In the BCS limit (Sec- 
tion VI) we will find that the quasiparticle scattering con- 
tribution gives the dominant contribution, while in the 
BEC limit (Section Vll) it is the zero point motion of 
the collective modes. The role of the convergence factors 
is explained in more detail in Appendix B, and further 
insight will also be found in the BEC limit. 



V. RESULTS FROM MEAN FIELD THEORY 
PLUS GAUSSIAN FLUCTUATIONS 

Once the thermodynamic potential is obtained, we can 
find the chemical potential as well as all thcrmodynam- 
ical variables of the system. We must, however, first 
determine the uniform, static gap parameter Aq. From 
eq. 124]) . we see that the Aq used in the expansion is 
the one that, for a given chemical potential, satisfies the 
mean-field saddle point equation ()lip , around which the 
action S is expanded to quadratic order. Thus, the gap 
and number equations 

SSo/SAo = and n = -dQ/dfi (28) 



FIG. 3: Ground state energy per particle £ in units of the 
non-interacting result 3e//5 as a function of — l/(fcFas). The 
difl^erence between the mean field result (dashed line) and the 
Gaussian fluctuation calculation (solid line) is small and more 
clearly shown in Fig. [6] (BCS limit). Fig. |1] (BEC limit) and 
Table I (unitarity). 

constitute the simplest theory which goes beyond the 
mean field approach and is consistent with Goldstone's 
theorem (see below). As we shall see in this and the 
next three Sections, this approach leads to very useful 
results and insights. We note that even though the sad- 
dle point gap equation pT|) used here retains its mean- 
field form, the values of Aq and ^ obtained from the 
simultaneous solution of (^5]) will deviate significantly 
from the mean field results (which are obtained using 
r^o of eq. p6|) in the number equation). Moreover, as 
we show in Appendix G, an identical approach leads to 
the known results in a different problem, that of a di- 
lute repulsive Bosc gas. In Section IX, we will analyze a 
dilfcrcnt scheme with a modified gap equation which in- 
corporates the self-consistent feedback of Gaussian fluc- 
tuations in the calculation of the saddle point, and show 
that it fails in some important aspects as an appropriate 
theory throughout the crossover. 

In this Section we present the results on the follow- 
ing quantities across the BCS-BEC crossover obtained 
by adding Gaussian corrections to mean field theory: (i) 
the gap parameter Aq (ii) the chemical potential /i (iii) 
the ground state energy £, (iv) the speed of sound c^, 
and (v) the Landau critical velocity. In the next three 
Sections, we will discuss the asymptotic results in the 
BCS and BEC limits and detailed numerical results at 
unitarity. 

In order to obtain Aq and fi from we solve the 

gap equation (fT4|) for Ao(^) together with the number 
equation written as 

Note that the thermodynamic /i-derivative (keeping vol- 
ume and r = fixed) in must take into account the 



/i-dcpendcnce of the saddle point Ao(yu). (An analogous 
point for the possibly more familiar case of the dilute 
Bose gas is emphasized in Appendix G). 

To solve the above number equations we must nu- 
merically evaluate ilg[fi, Ao{fi)]. Even though the real- 
frequency representation (j27p gives physical insight, we 
find it simpler to numerically evaluate fig on the Mat- 
subara axis, as described in Appendix C. Finally, we cal- 
culate F{p) = Qq + fic,[/i, Ao(/i)] + fin and look for an 
extremum (maximum) as a function of the chemical po- 
tential /i. 

The gap Ag and the chemical potential fi are plotted as 
a function of —{kpas)^^ in Fig. [21 where the dashed line 
is the MF value for comparison. As expected, the inclu- 
sion of fluctuations reduces the value of Aq. We note that 
the auxiliary field Aq continues to determine the energy 
gap Eg = Ao for > and Eg = ^l^p + Aq for /i < 



(just as in MF theory), so long as we ignore the feedback 
of the fluctuations on the single-particle propagator. 

The ground state energy of the system is obtained from 
the thermodynamical potential using £ = ri(T = 0) -I- fin 
is plotted in Fig. [3l We see that although the differ- 
ence between the MF and Gaussian results is quite small, 
fluctuations reduce the ground state energy through the 
entire crossover. The quantitative superiority and the 
physical insights of the Gaussian results are discussed in 
detail later: see Fig. [6] for the BCS hmit, Fig. [8] for the 
BEG limit and Table I for the results at unitarity. 

We next compute the speed of sound through the BGS- 
BEG crossover. First, we emphasize that Goldstonc's 
theorem is necessarily obeyed by the theory defined by 
(PS)) : this is in contrast to the self-consistent calculation 
to be described in Section IX. The existence of a zero 
energy Goldstone mode is guaranteed by the form of the 
gap equation (fTT|) which implies that DetA/ (q = 0, a; 
0) = 0. To see this fact, note that we can write 



DetM(0, 0) = 



^DetG 



(30) 



and the saddle-point condition is 1/g = — J2k DetG. 

The collective mode spectrum has the form ujo{q) = 
Cs|q| for g where Cg is the speed of sound. We cal- 
culate Cs following the approach of we include the 
expressions here to correct a typographical error [33| in 
that reference. Rotating the frequency from the Mat- 
subara axis to the real line {iqi — > —tu) and expanding 
to quadratic order in both momentum and frequency, 
we get Mii(-cj,q) = (A + 2Buj + (C + Q)|q|2 - {D + 
i?)u;2)/2, M22(-c^,q) = A/ii(w,q) and Mi2(-Cc;,q) = 
(A+ (C- Q)|q|2 - {D-R)uj^)/2. Here A = Ek ^o/^E^ 

B = Eke/2£;^ c - Ek{(i - m^h^ - [i - 10^(1 - 

X)]Y}/^E\ D ^ Ek(l - X)/%E^ Q ^ EkiCM - (1 - 
3X)Y}/8E^, and finally R = Ek with the nota- 




FIG. 4: The Landau critical velocity as a function 
of —l/{kFas) in the BCS-BEC crossover is given by 
min{£ {k) / k} . The black lines represent the speed of sound 
Cs, and the gray lines the pair breaking estimate obtained 
using (|32p . In each case the dashed line is the result using 
the MF gap and chemical potential, while the solid line is 
the result obtained after inclusion of Gaussian fluctuations. 
The black data points (with error bars) are the experimental 
results from ref. [s^. The open squares are the result of a 
Bogoliubov-deGennes vortex calculation from ref. [2^ . 



tion X = Al/E^, and Y = [kp/Sm^. We thus obtain 

Cs = VQ/[ByA + R]. (31) 

The results for the speed of sound across the BGS-BEG 
crossover are shown as the black curve in Fig. |4l The 
solid line is the result obtained after inclusion of Gaussian 
fiuctuations, while the dashed line is the result using the 
MF values for Ag and fi. The other curves shown in this 
Figure are discussed below. 

To conclude this Section, we turn to the calcula- 
tion of the Landau critical velocity Vc as a function of 
— l/(fcFas)- As we have seen, most observables - gap, 
chemical potential, ground state energy, speed of sound - 
are monotonic functions oil/kpas through the crossover. 
The same is true of the transition temperature Tc which 
shows a slight maximum near unitarity [isj but is essen- 
tially the same, of order 0.2e/ for all positive scattering 
lengths, i.e., to the BEG-side of unitarity. In other words 
there seems to be nothing particularly dramatic about 
the properties of the most strongly interacting unitary 
regime. However, as first pointed out in ref. [28| based 
on the study of the current flow around a vortex, the 
critical velocity as a function of l/kpas has a strongly 
nonmonotic behavior through the crossover with a pro- 
nounced peak at (or close to) unitarity. The reason for 
this behavior is that very different excitations are respon- 
sible for the destruction of superffuidity [2^, [2^: breaking 
of pairs on the BGS side and generation of phonons on 
the BEG side of unitarity. 

To understand this better and to compare with recent 
experimental data, we use the Landau criterion which 
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gives an upper bound on the critical velocity of the form 
Vc = min{£'(fc)/A;} where £{k) is the energy of an ex- 
citation carrying momentum k. We separately consider 
single-particle (fermionic) and collective (bosonic) exci- 
tations. For single-particle (s.p.) excitations, the excita- 
tion energy is -Bk = \/(ek - A*)^ 
to the Vr estimate 



Aq which then leads 



s.p. 



1/2 



(32) 



As wc shall see, this is most relevant on the BCS side of 
unitarity, and in the BCS limit fi 3> Ao, it simply reduces 
to the well-known result for pair-breaking Vc — Ao/kp. 
The pair breaking estimate of eq. ((32)l is plotted in gray 
in Fig.m The dashed gray line is the result using the MF 
gap and chemical potential, while the solid gray line is the 
result obtained after inclusion of Gaussian fluctuations. 

For collective excitations, we find that the Landau crit- 
ical velocity is given by the slope of the tangent to the 
LUQ{k) curve. Since there is no roton dip for the super- 
fiuid Fermi gas (as already remarked), one simply gets 
the speed of sound: 



(«c) 



coll. 



(33) 



The actual critical velocity is then bounded above by 
the minimum of the single particle {vc)^ p and collective 
ivc)^Q\i -An estimate of the Landau critical velocity at the 
mean field level, which corresponds to the dashed curves 
in Fig. 4, was given in [2^. We find that quantum fluc- 
tuations lead an appreciable reduction in Vc, as seen in 
the full curves in the figure. We also plot the results of 
a recent experimental study of the critical velocity [sO], 
for which Vf = 30 mm/s. The theoretical predictions of 
the nonmonotocity of Vc with a peak around unitarity 
and their experimental confirmation show that the uni- 
tary Fermi gas is the most robust superfluid in the entire 
crossover. 



VI. BCS LIMIT: HARTREE SHIFT AND FERMI 
LIQUID CORRECTIONS 

We now describe in detail the BCS limit solution 
for 1/kpas —oo. We will show that the collective 
mode contribution in cq. (|27p is very small because of 
phase space restrictions, and the dominant correction to 
fl{T = 0) comes from virtual scattering of quasiparticles. 
We find that the Gaussian theory recovers the well-known 
"normal state" correction to the ground state energy of a 
dilute Fermi gas originally studied by Huang, Yang and 
Lee and by Galitskii. The leading term in this correc- 
tion, which is of order kpas, is the Hartrce shift of the 
ground state energy. We note that the BCS mean field 




FIG. 5: Spectrum of excitations which contribute to the Gaus- 
sian correction to the thermodynamic potential (|27|) . The re- 
sults shown correspond to the BCS regime l/{kFas) = —2, 
with ^ = 0.867e/ and Aq = 0.0311 e/. The full line shows 
the collective mode dispersion u!o{q) (pole of 1/DetM) and 
the shaded region is the two-particle continuum (branch cut) . 
The dashed lines uj22{q) and Wii(g) are the zeros of , M22, and 
Mil respectively, which coincide in the particle-hole symmet- 
ric BCS limit. The inset shows long- wavelength, low-energy 
spectra. 



ground state energy differs from the free Fermi gas by 
a condensation energy of order Aq/ej and represents an 
exponentially small correction of order cxp{~l/ kp\as\) 
relative to e/. In contrast, the Gaussian fluctuation con- 
tributions will be found to be power-law corrections in 
kplasl- 

The BCS limit is characterized by an exponentially 
small gap Aq and Cf. The spectrum of collective 
excitations found in the BCS regime is shown in Fig. [5]for 
l/ikpas) = -2 with n = 0.867e/ and Ao = 0.0311 e/. 
In this case, at 5 = the continuum starts at a frequency 
equal to 2Ao, and the collective modes arc restricted to 
a small frequency interval, magnified in the inset. Due 
to particle-hole symmetry in the BCS limit, the zeros 
LO22 and u)ii coincide and hence do not contribute to the 
energy (j27p . The speed of sound in this limit becomes 
Cs — Vf /VS [sHi as we show at the end of this section. 

The MF ground state energy of the superfluid state 
is given by the well known BCS result £q = 3nef/5 — 
(3nAo/8e/). We first show that the contribution of the 
zero-point motion of the collective modes to f2 is expo- 
nentially smaller than the (already small) MF conden- 
sation energy, and may be neglected. The momentum 
Qc where the pole hits the continuum at 2Ao is given 
by Qc ~ Ag/w/ ~ where ^ is the correlation length. 
Thus the phase space available for the collective mode 
contribution is tiny because qc ^ kp. The contribution 
per particle to (^7)) coming from the poles is seen to be 
~ f-fikpS)"^ ^ e/(Ao/e/)^ which is negligible. 

We next turn to the continuum contribution to (|27|) . 
In the BCS limit it is justified to set Ag = here to 
get the leading order terms in the ground state energy. 



Any corrections due to non-zero Aq are at least down 
by a factor of 0(Ag/e/) and thus negligible. As already 
remarked, the ground state is a superfluid which leads to 
a MF energy reduction of ©(Aq/ej), with respect to the 
energy that we will calculate. This exponentially small 
contribution is vital to get a stable ground state, but once 
this is done, we can set the gap to zero in computing 
leading order corrections to the ground state energy, as 
explained in Appendix D. 

In the BCS limit, our result for fig is exactly of the 
form of the well-known results for the dilute Fermi gas 
in its normal state f26j . but with as < 0, as shown in 
Appendix D. We find that the total energy is given by 



nef 



3A| 

8.^ 



(34) 



-—kpas 
37r 



4(11 - 21n2) 



357r2 



{kpasf 



and the chemical potential 
^ = 1 - 0(Ao/6/)2 



(35) 



4 , 4(11 -21n2),, ,2 
+ ir^Fas + , -{kFasY + ■ • ■ 

37r lOTT^ 

The first order term is the Hartree term while the second 
order contributions have the same form as those obtained 
by Huang, Yang and Lee and by Galitskii 24 , 2^ for the 



dilute Fermi gas, except that in our case 05 < 0. 

It is worth commenting that, although it is customary 
to think of the Hartree term as a "mean field shift" , it 
arises in our approach as the first term in the fluctua- 
tion correction to the saddle point thermodynamic po- 
tential in the BCS limit. This is also seen clearly from 
Fig. [To] in Appendix D where the first diagram is clearly 
the Hartree term. Once we move away from the BCS 
limit toward unitarity, the gaussian fluctuation contri- 
bution is still well defined even in absence of a small 
parameter, however a "Hartree term" becomes hard to 
identify. 

In Fig. [6] we plot the energy per particle as well as the 
chemical potential as function of interaction in the BCS 
regime: —2 < (fc^as)"^ < —1. We have also included 
the mean field result and the asymptotic values given 
by ([M)) and ([55)) . As we can sec, quantum fluctuations 
reduce the energy and the chemical potential, in a way 
that is consistent with the corrections obtained from our 
analysis. 

In order to calculate the speed of sound, one needs 
to look at the slope of the pole dispersions, which in a 
superfluid is given by c = \J psj where ps is the su- 
perfluid density and k = mdn/dp is the compressibility 
of the system. In a homogeneous system, Galilean in- 
variance implies that ps = 71 at T = 0. Keeping only 
the Hartree term in the formula for the chemical poten- 
tial, we obtain to linear order in the scattering length 
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FIG. 6: Energy and chemical potential in the BCS regime, for 
-~{kFas)~^ between 1 and 2. The solid lines are our calcu- 
lations including Gaussian fluctuations, the dashed lines the 
mean field values, and the dotted lines the results of formulae 
(|34p and (|35p . in the top and bottom panels, respectively. 



Cg ~ w//V3(l -|- ^kpas), which was first obtained for 
BCS superconductors by Anderson (35j . 

We conclude this Section with a comment on the 
Gorkov Melik-Barkhudarov correction [s^ which enters 
the pre-cxponential factor in the BCS expression for the 
gap. It arises from the renormalization of the effective at- 
traction by particle-hole excitations in the medium. Such 
effects are not accounted for in our theory which effec- 
tively considers only particle-particle channel diagrams. 
At the present time we do not know of any theory which 
recovers this correction in the extreme BCS limit and 
shows how this correction evolves through unitarity. 



VII. BEC LIMIT: DIMER SCATTERING AND 
LEE- YANG CORRECTIONS 

We next describe in detail the BEC limit solution for 
1 /kpas +00. We will see that the zero-point motion of 
collective modes entirely dominates over the continuum 
contribution in the thermodynamic potential. We will 
find that this leads to two important effects in the BEC 
limit: first, a reduction of the effective dimer-dimer scat- 
tering length relative to its mean field value of 2as and 
second, a Lee- Yang correction to the equation of state 
of the dilute gas of dimers. While our theory is able to 
obtain both these effects semi-quantitatively, it does not 
give the exact asymptotic answers. The scattering length 
for bosonic molecules, or dimers, is found to be ~ 0.55as, 




FIG. 7: Spectrum of excitations which contribute to the Gaus- 
sian correction to the thermodynamic potential (|27p . The re- 
sults shown correspond to the BEG regime with l/{kFas) — 
2, with n = -3.967eF and Ao = 1.025eF.. The fuU line 
shows the collective mode dispersion u)o{q) (pole of 1/DetM) 
and the shaded region is the two-particle continuum (branch 
cut). The dashed line iJ22{q) are the zeros of M22. We do 
not show ujii{q), the zeros of Mn, because they coincide with 
the continuum in the BEG limit. The inset shows the long 
wavelength, low energy spectra. 



while the exact solution of the four-body problem yields 
0.605 [S^l, and the coefficient of the Lee- Yang correction 
is only 6% smaller than the exact result. 

In the BEC limit the chemical potential is large and 
negative and, to leading order, goes to one half of the 
binding energy of the molecules: /x = —E{,/2 where Eb = 
l/mog. The spectrum for collective excitations is shown 
in Fig. [7] for 1/kpas = 2. The two-particle continuum 
then sits at a very high energy; at |q| = it begins at 
an energy of 2y^/xP + Aj « Ef,. The virtual scattering 
of the very high energy fermionic excitations makes a 
negligible contribution to the thermodynamic potential 
in the molecular BEC limit. 

The Gaussian contribution is then entirely dominated 
by the low frequency collective mode uJo{q), which is the 
Bogoliubov excitation of the molecular Bose gas, and 
the mode 1022(1) coming from the convergence factor. 
These modes are shown in more detail in the inset of 
Fig. [71 Note that for ^ < 0, wn, the zeros of Afn, 
are pushed to the continuum and do not enter the cal- 
culation. The thermodynamic potential then simplifies 
to Qg ~ l/2X]qN(q) - t^22(q)]- In the BEC limit the 
pole always remains below the continuum for all q and 
the phase space for the zero-point oscillations extend for- 
mally to g = cx). 

This raises the question: how can we get a finite an- 
swer for the sum in fig? The dispersion of the pole is 
given by the standard Bogoliubov expression wo(ct) = 
■\/c^q^ -f- (q^/2m;,)^, where = 2m is the mass of the 




l/(kpa) 



FIG. 8: Energy and chemical potential in the BEG regime, for 
(kpas)'^ between 1 and 2. The solid lines are our calculations 
including Gaussian fluctuations, the dashed lines the mean 
field result, and the dotted lines the leading order result for 
a dilute Bose gas with an effective repulsion ai, = O.Gas- 



of this dispersion is just the kinetic energy of the bosonic 
molecule, which should not be part of the zero-point mo- 
tion of the fluctuations. This is exactly where the conver- 
gence factors come in. One finds that (jJ22(q) q^/2m(, 
for large q and cancels the contribution of the free bo- 
son dispersion. Thus the convergence factor gives rise to 
manifestly finite results by eliminating the free particle 
dispersion of the pole spectrum from contributing to the 
zero-point motion of the phase fluctuations. (See Ap- 
pendix G for an analogous discussion for the dilute Bose 
gas.) 

We now calculate the leading order correction to the 
mean field thermodynamic potential in the BEC limit, 
which will in turn determine the effective interaction be- 
tween the bosonic molecules. The idea is to use an ex- 
pansion in the small parameter Ao/|/i| ^ 1 in the BEC 
limit. The calculation is most easily done on the Matsub- 
ara axis as detailed in Appendix E. Here we only quote 
the leading order result: 



2567r 



(2m) 



3/2 



A* 
^0 



3/2 



(36) 



bosonic molecule as shown in ref. 1J|. The large q limit 



where we have included the factor of 2567r in order to sim- 
plify later expressions. The m, Aq and /i dependence can 
be determined analytically and the dimensionless prefac- 
tor a = 2.61 has to be evaluated by a numerical integra- 
tion, as shown in Appendix E. 

To find the effective scattering length between the 
molecules we calculate the shift in the bosonic chemi- 
cal potential as a result of interactions. Toward this end 
we proceed as follows. We expand the gap equation ^ 



11 



in powers of Ao/|^| <C 1 to recover the Gross-Pitacvskii 
equation for bosons. We then get 



-= 



as 



(37) 



As noted earher, the leading order result is /i = 
— l/(2ma|). The above result allows us to relate the 
next order correction in fi to Aq. We find 



-l/(2mas) + Sfi, 



Sfj. = magAQ/4. 



(38) 



Next we determine Aq from the number equation (|29p . 
We find: 



n = umf 



dria 



9A| 



(39) 



where timf = Alm?/"^ / {A^/2■K\^l\^/'^) is the MF result 
in the BEG limit. (See Appendix G for an analogous 
discussion of keeping track of the /i dependence of the 
saddle point for a dilute Bose gas.) Using ([36|) we get 
dVlg/dn = (a/2567r)(2m)3/2(3/2)A;^/|^|5/2 and the last 
term in ^ is given by {a/?,2Ti){2'mf/'^Al/{\^L\^/'^maj^). 
It is easily seen that —dQ.g/dii ~ {kpasY ^-^^d can be 
neglected in the BEG limit. Now, using = l/(2ma|.) 
to leading order, we get 



A^ 

^0 



1 



kpas \1 



(40) 



This leads to (5/x = (2/37r)e/A:Fas(l/(l + «)). Gomparing 
the chemical potential for the bosons fXb — 25/x with the 
weakly interacting Bose gas result /i& = AirabUb /mt with 
Ub = n/2 and m?, = 2m, we get the effective scattering 
length for the bosons to be 



Of, = 2as/(l + a) ~ O.SSag, 



(41) 



using our numerical result a = 2.61. This result for Oh 
is identical with the one obtained by Hu et al. [l^l; see 
Section X for further discussion. 

Going beyond the leading order term we find that the 
next order correction to the chemical potential is of order 
{nbaiy/'^{nbab/mb), which has the same form as the Lee- 
Yang corrections for a weakly repulsive Bose gas [l^l . To 
see this we analyze our numerical results for O in the 
BEG limit as follows. Scaling out energies with l/(2ma|) 
and lengths with as, we find that we can fit to the 
functional form Og = {l/2mal)a^^[A{Sfi)^ + B{Sii)^^^ + 
. . .] where = 2mag^ — 1. Solving for the molecular 
chemical potential /Jb ~ 26fi we obtain 



IJ-b 



Annbab 
rub 



1 + 7 



32 
3V^ 



(nbal)'^' 



(42) 



where we find that the coefhcient 7 = 0.94 is 6% smaller 
than the Lee- Yang result 7=1 (see Appendix G). 



In Fig. [8] we show the total energy and the chemical po- 
tential for interactions in the BEG regime 1 < l/kpas < 
2. For comparison we also show the mean field results 
(dashed lines) and the leading order result for a gas of 
bosons interacting with ab = 0.6as (dotted lines). 



VIII. UNITARITY 

At unitarity {1/kpas = 0), there is no small parameter 
and the problem can only be solved numerically. At this 
point both the pole and the continuum corrections are of 
comparable magnitude in the thermodynamic potential. 
We present the results for the gap, chemical potential, 
ground state energy and the speed of sound and compare 
with quantum Monte Garlo and experimental results in 
Table D 

At unitarity \as\ diverges and this leads to the con- 
cept of universality, i.e., all energies scale with the Fermi 
energy e/ and all lengths scale with kp^. A conse- 
quence of universality is the relation between the ground 
state energy per particle and the chemical potential /i = 
{5/3){£ /N), which acts as a check on our numerical cal- 
culation. The ground state energy is generally written in 
terms of the non-interacting energy £/N = (1 -I- /3)3e//5. 
(Note that in this Section /3 is used to denote the uni- 
versal interaction correction to the ground state energy, 
and not the inverse temperature). We obtain a numeri- 
cal value oi (3 = —0.598(1). The mean field theory gives 
f3 = —0.41 [3], while quantum Monte Garlo methods 
give a /3 of —0.56 The experimentally obtained val- 
ues range from -0.68 to -0.49 0, 0, [l3|- We thus see 
that at unitarity the Gaussian quantum fluctuations - 
zero-point motion of collective modes and virtual quasi- 
particle scattering - account for most of the difference be- 
tween the exact ground state energy (i.e., that obtained 
from QMG or experiments) and the simple mean field 
estimate. 

The speed of sound is obtained from the dispersion of 
the pole Lo = Cs |q| at small momenta, or alternatively in 
our theory from eq. (PT|) . However, we can also calculate 
the speed of sound once we know the equation of state. 
Using that d^L/dn = (2/i)/(3n), we arrive at the expres- 
sion Cs/vj = ^(1 + /3)/3. Our theory predicts a speed 
of sound at unitarity of Cg = 0.37 w/, using either one 
of the mentioned methods. For comparison, the answer 
14| obtained by using the mean field gap and chemical 
potential is OAAvf. The Quantum Monte Garlo estimate 
is Cs = 0.38v/ while the experimentally measured value 
of the speed of sound at unitarity is 0.38u/ [STf . 

We also show in the last column of Table I various 
gap estimates. The inclusion of fluctuations reduces the 
value of Aq relative to the MF estimate. As already 
noted above, Aq continues to determine the energy gap 



Eg ~ Aq for /i > 0, so long as we ignore the feedback of 
the fluctuations on the single-particle propagator. We do 
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not include an experimental value for the energy gap as 
we believe that it is not clear how to quantitatively ex- 
tract this from rf spectroscopy data, taking into account 
the interactions between atoms in the final and initial 
states 
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IX. SELF-CONSISTENT FEEDBACK OF 
GAUSSIAN FLUCTUATIONS ON THE 
SADDLE-POINT 

In this Section we describe how one can include the 
feedback of the Gaussian fluctuations on the saddle-point 
equation in a self-consistent manner. We show here that 
it is straightforward to accomplish this using a fluctu- 
ation formalism similar to the one used in Section IV. 
However, one finds, quite generally, that Goldstonc's the- 
orem is violated if one uses such a Cartesian represen- 
tation for the fluctuations, as soon as one modifies the 
saddle point equation. We next show that a polar rep- 
resentation of the fluctuations in terms of the amplitude 
and phase of the auxiliary field allows one to recover the 
Goldstone mode, even when the saddle point condition is 
modified away from its mean field form. There is still a 
problem with obtaining an ultraviolet convergent expres- 
sion for the thermodynamic potential in terms of the fluc- 
tuations. The polar representation respects Goldstone's 
theorem at low energies but has unacceptable high energy 
properties, while the Cartesian representation violates 
Goldstone in the infrared but is well controlled in the 
ultraviolet. We thus construct a hybrid representation 
which interpolates between the polar in the infrared and 
the Cartesian in the ultraviolet, and compute the ther- 
modynamic potential. In the end, we are not convinced 
from the solutions of the new gap and number equations 
that the self-consistent theory is worth the effort. In fact 
we find results which are not an improvement relative 
to those presented in Section VIII at unitarity and the 
theory has problems in the BEG limit. 

Formally the calculation proceeds in much the same 
way as in Section IV with one important difference. We 
again write A (a;) = Aq + ri(x) where Aq is a real num- 
ber, which is the (x, T)-indepcndcnt part of A(a;), and 
r]{x) are the complex fluctuations about it. We call this 
the Cartesian representation of fluctuations to be con- 
trasted with the polar representation to be introduced 
below. The main difference with Section IV is this: Here 
Aq does not follow the mean field gap equation and its 
value will be determined only after including the effect of 
fluctuations, as explained in detail below. In this sense, 
the ?7's are not the fluctuations about the mean field sad- 
dle point, but rather about a uniform static value Ag 
which will be determined self-consistently after integrat- 
ing out the fluctuations. 

We again find that to order rj^ we get the action 



((12)) even though Aq is not set to its mean field value. 
The Gaussian piece too has the same form as ([20]) . We 
emphasize that there is no linear term in rj in eq. (jl9p . 
despite the fact that we are not expanding around a sad- 
dle point. The reason for the absence of a linear term 
is that such a term would be proportional to rj^q = 0). 
However, r]{q = 0) = since the uniform (q = 0), static 
(iqi = 0) piece of of A{q) is described by the (as yet 
undetermined) Aq. 

Integrating out the Gaussian fluctuations, we obtain 

Z ~ j '^^o J DrjDrj^ exp (-5o - Sg) (43) 

= y"dAoexp[-S'eff(Ao)] (44) 

with the effective action 

^cff = ^0 + (1/2) lnDctM(<7). (45) 



where was defined in eq. p2|) |39j . Using the con- 
vergence factors described in Section IV and Appendix 
B we get the final result 

^^f 4E-lnG-(.) 
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(46) 



Sa — So + Sq 



where Sq has the mean-field-like form 



The parameters Aq and fi arc then fixed by solving the 
gap equation given by 

SS,a/S Ao - SSo/SAo + SSg/6Ao = (47) 

and the number equation given by 

n =^ -dn/dfi = -^[So/(3]/^^i - d[Sg/f3]/dn (48) 

The theory developed above has a serious problem: 
there is no zero energy Goldstone mode in the system. 
To see this recall cq. (PH)) for DetM(0, 0) (which contin- 
ues to be valid here). Wc immediately see DctA/(q = 
0, w = 0) 7^ unless l/.g = - X^kDctG. The last condi- 
tion is however the mean field gap equation which is not 
satisfied by solutions Aq and ^ of the self consistent gap 
equation (|i71) and the number equation We thus 

see that in the Cartesian representation, the Goldstone 
mode is lost as soon as one moves away from the mean 
field saddle point. 



Amplitude and Phase Fluctuations in the 
Self-Consistent Theory 

How can we restore the gapless Goldstone mode in 
a self-consistent Gaussian calculation? This can be 
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achieved by using a polar representation for the fluctua- 
tions in terms of amphtude and phase: 



A{x) ^ Ao[l + A(a;)]e 



ie{x) 



(49) 



in place of the Cartesian representation A(x) = Ao+rj{x) 
used above. We wiU first show that the phase excita- 
tions 6{x) are necessarily gapless in the long-wavelength 
limit, even when the saddle point shifts away from the 
mean-field value. However, we will find that there is a 
price to pay for obtaining the correct low-energy, small- 
|q| physics. The high-energy, large-|q| behavior of the 
amplitude-phase fluctuation propagator has unphysical 
properties, and finally we will be forced to an interpo- 
lation scheme between the polar representation at low 
energies and the Cartesian representation at high ener- 
gies. 

Working with the amplitude X{x) and phase 9(x) we 
get 



Z = J dAo J D\D9J cxp{-Sao,x,o)- 



(50) 



A detailed derivation of the results stated here is given 
in Appendix F. As shown there, the action is the sum of 
two terms 



Sao.xm — So + Sg, 



(51) 



where 5*0 defined in eq. ([T2j) has the mean-field form, and 
the Jacobian J of the transformation is approximated by 



j=nA 



(52) 



This is the same approximation used in Section IV A of 
[40| . We will see later that the contribution from J is 
exactly canceled by another contribution [see below eq. 

The Gaussian term Sg can be written as 



We use the notation Sg here to distinguish it from the 
Gaussian action Sg in the Cartesian case ()20p . The in- 
verse fluctuation propagator D is given by 

^0 , ^0 



Dii 
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J2 trGo(fc)criGo(fc-fg)CTi 



D22 = |- V trGo(fc)a3 



-trGo(fc) (igjCTa - S^)Go{k + q){iqiCT3 - S^) 



D12 = ^ 5] tYGoik) {tqi<j3 - SOGoik + q)ai 

D21 = D12 



(54) 



where the Pauli matrices (Ji operate in Nambu space and 
= Ck+q — = ek-(-q — Ck- 

Infrared behavior: The long-wavelength, low-energy 
limit of the amplitude and phase fluctuations described 
by eq. ((54|) have the following properties, (i) Di2(q = 
0, iqi = 0) = 0, so that the amplitude and phase modes 
decouple in the q ~ limit, (ii) The 'D22{q)\d{q)\'^ term, 
upon transforming to space-time, has form /9s|V0p/2 — 
K{d9/dt)^ + . . ., where ps is the superfluid density and k 
the compressibility. In particular, we note that D22(q = 
0, iqi = 0) = for any choice of Aq and /i, so that the 
phase mode is gapless in the long wavelength limit. Thus 
Goldstone's theorem is respected even when one moves 
away from from the mean field saddle point, in marked 
contrast to the case of Cartesian fiuctuations. 

Now, it would seem that fiuctuations in the amplitude- 
phase representation appear to solve all our problems. It 
is tempting to argue that one can simply integrate out the 
A and 9 fields in eq. ([SO)) and obtain an effective action 
which is the analog of cq. (|45)) with J2qiqi InDetM(g) 
replaced by J2q iq, lnDctD((7). However the situation is 
not so simple. As we show next, the high-energy be- 
havior of DetD((7) is such that the required Matsubara 
sum diverges, and there is no analog of the convergence 
factors in eq. ((46|) . 

Ultraviolet behavior: Wc find it useful to rewrite 
the D-matrix in a form which permits us to better un- 
derstand its high energy properties and also to see its 
relationship to the M matrix used to describe the fiuc- 
tuations in the Cartesian representation. Omitting the 
rather lengthy algebra involved (which is sketched in Ap- 
pendix F), we find that 



Dii 



Do 



2 ^ 

^ k,iA:„ 



fz-io /-lO/ 1/^0 (-'0/ _j_ ono con 



^0 



[^2265^ -I- G22G51 — 2G52G"2 — 2DetG°] , 



k,2/c^ 



D12 = -D21 = ^ E [G22G!'i - G^liG%] , (55) 

k,ife„ 

where we have used the notation G° = G°(fc) and G°' = 
G''(A;-|-g). Note that both the properties discussed below 
eq. ([54)1 - the decoupling of the amplitude and phase 
modes at g = and the Goldstone theorem - are also 
evident in the new expression for the D matrix. 

From the expansion of the order parameter field in 
eq. ip^l to linear order, we see that the fluctuations of the 



order parameter are of the form A = AqA and 9 = A06'. 
It is useful to rescale the fluctuation fields to A and 0, 
so that D = D/Aq. This leads to a factor of 21nAo in 
the action which exactly cancels the factor J in Eq. ([5^ . 
In this form it is also easier to make connection with the 
Cartesian 77 fields, as shown below. 

At sufficiently high energy and/or short distance 
scales, the system must "look normal" (i.e., non- 
superfluid) and the natural variables to describe the fluc- 
tuations are the Cartesian 77's: 

In this basis the matrix D ^ L = WDW^ where 
Lii = i(D22 + Dii - 2iDi2) = Mil - ^ 

1 ~ ~ ~ X 

L22 = 2 (D22 + Dii + 2iDi2) = M22 - Y 



2 




q/kp 

FIG. 9: Spectra used in the self-consistent calculation of 
Gaussian fluctuations. As explained in the text, the method 
requires to modify the contribution of the poles of the fluc- 
tuation matrix at low energies, from that of the cartesian 
representation to that of the amphtude-phase representation. 
The interpolated values used in the calculation are included 
in the dashed line. 



L12 =L21 = -(D11-D22) =Mi2 + y (57) 

where the M matrix was defined in eqs. (|2f I22p and 
X = l/9 + Y. DetGo(fc) = l/g-Yl 1/{2E^)- (58) 

k k 

We now see that, insofar as fluctuations about the mean 
field saddle point are concerned, the L and M matrices 
are identical. This follows from the fact that X = Q when 
the saddle point equation has the mean field form. (This 
was the case in the calculation described in Sections V 
through VIII, even though Ao and ^ did not have their 
mean- field values.) 

Conversely, if we look at fluctuations about a saddle 
point defined by an equation which does not have the 
mean field form - which is the case here - then X ^ {) 
and the inverse fluctuation propagators M (directly ob- 
tained in the Cartesian representation) and L (obtained 
by transformation from the polar D to the Cartesian rep- 
resentation) necessarily differ. It is only L, derived from a 
polar representation, that respects Goldstone's theorem. 
The presence of the X factors which ensure the Goldstone 
mode in the infrared, however, spoils the ultraviolet be- 
havior of L and prevents one from getting a convergent 
answer for X^iq, lnDetL(q, iqi). The mathematical analy- 
sis showing this difficulty is sketched in Appendix E; here 
we give a simple argument which indicates the problem. 

After analytic continuation from iqi u) + iO"*" we find 
that in the w — 00 limit, the L matrix looks schemat- 
ically like L22(q, ljj) ^ —l/as + i-\/M ~ + . . . and 
Lii(q,w) -1/as -f- yj\uj\- X/2 + ... and Li2(q,cj) 
|i^|-3/2 _|_^y2. (We omit multiplicative constants here in 



various terms and simply focus on their dependence on 
as,uj and X). The presence of the X/2 factor in the L12 
fundamentally changes its asymptotic behavior from the 
jo;!^'^/^ in the M matrix to a constant X/2. This leads 
to convergence problems discussed in Appendix F. 

Results of the Self-Consistent Calculation 

We resolve the problem described above by choosing 
a scheme that interpolates between the polar represen- 
tation in the infrared and the Cartesian representation 
in the ultraviolet. We define an "interpolating" collec- 
tive mode energy w^(q) = wo(q) + /(q)(Lj^(q) - wo(q)) 
where we choose /(|q|) = 1 — cxp (1 — gc/|q|) for |q| < Qc 
and /(|q|) = for |q| > qc, with qc ~ kp (the only scale 
at unitarity). This formula goes smoothly from Wq ~ Wq , 
the polar result at small g to Wq ~ wqj the Cartesian re- 
sult at large g, as shown in Fig. [9l Since the continuum 
contributions occur at high enough energies (at unitar- 
ity) we have left these (Cartesian representation) con- 
tributions untouched. Operationally, we implement this 
by adding the following term into the thermodynamical 
potential 

&^sc = i^/(q) [^0(9) -'^o(q) -'^^2(q) +'^22(q) 

q 

+c^fi(q)-c^ii(q)]. (59) 

Using this method and solving the modified gap and 
number equations (|47l |48|) we obtained that at unitar- 
ity the chemical potential is /isc = 0.35 e/ with a gap of 
Ao,sc = 0.68 ey. These values compare rather unfavor- 
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ably with the quantum Monte Carlo values as well as the 
experimentally measured values. 

We next show that the self-consistent calculation has 
serious problems in the BEC limit: the effective interac- 
tion between the bosons is attractive and the system is 
thermodynamically unstable! Clearly this is an artifact 
of the modified gap equation. In parallel with the anal- 
ysis in Section VII we can show that, just as in eq. ([36]), 



ng ~ -a(2m)3/2A;^/|^|3/2/(2567r). Although the values 
of Aq and will change because now we are using a 
modified gap equation, the dimensionless constant a re- 
mains the same as before. As shown in Appendix E, it 
is given by a = 2.61. From the new gap equation wc can 
show that cq. ((37|) is now modified to 



— = V2m|/.j| 



as 



1 + (1 - Q') 



l6\^i\■■ 



from which wc find /i = — 1/ (2mag) + Sji with 



Sj-i = (1 — a)- 



-1.61 



(60) 



(61) 



Unlike the result ([38]) for Sfi in Section VII, we find that 
the self-consistent calculation yields (5/i < in the BEC 
limit. A reduction in ^ relative to the non- interacting bo- 
son value is equivalent to an effective attraction between 
the bosons or a negative compressibility. 



X. RELATION TO OTHER APPROACHES 

We now turn to a discussion of the relation of our work 
to that of other authors. First, the idea of writing the 
ground state energy density of a many-body system in 
terms of the zero-point motion of collective excitations 
(plasmons) goes back to early work on the electron gas 
using the "Random Phase Approximation" (RPA) 41 1. 
The RPA was generalized to the BCS supcrfluid in the 
early work of Anderson 35[, where the collective mode 
spectrum and its modification by long-range Coulomb in- 
teractions was discussed, but the question of the ground 
state energy density was not fully addressed as far as we 
can see. The inclusion of thermally populated collective 
excitations was central to the Nozieres-Schmitt-Rink the- 
ory of Tc in the BCS-BEC crossover 0, [H. In fact that 
was the dominant contribution on the BEC side of the 
crossover. The difference here is that we are looking at 
quantum corrections about the broken symmetry state 
where we have to deal with matrix propagators. 

Several recent works introduce a small parameter by 
hand; either by expanding in dimensionality around four 
or two dimensions [20| or by introducing a large number 
2A^ of fermion flavors with a Sp(2A^)-invariant Hamilto- 
nian [U . Our self-consistent calculation in Section 
IX is closely related to the "1/A^ expansion" approach. 
At zcroth order in 1 /N one obtains the mean field results 



and first order in 1/N gives the RPA or Gaussian fluc- 
tuations. The saddle point is then recalculated to low- 
est order in a 1/A^ expansion, with changes in the gap 
and chemical potential from their MF values obtained 
perturbatively in 1/iV, which is treated as a small pa- 
rameter. In practice the calculation is done to flrst or- 
der in 1/A^ (although in principle it could be done to 
higher order) and N is set equal to unity at the end. 
On the other hand, we keep = 1 throughout the self- 
consistent calculation. Thus the actual values of the Aq 
and /i obtained at unitarity, for instance, are quite differ- 
ent in our approach and in the large approach, even 
though if one was to set A^ = 1 throughout the equa- 
tions would look the same. One has to be rather careful 
about how various physical quantities are calculated in 
the 1/A^ expansion. For example, in the BEC limit we 
can show that = (47rn;,/mb) 2as(l — a/N) which is 
negative for A = 1 and would lead to a negative bulk 
modulus dfib/dub- On the other hand, the more natural 
quantity to compute in the grand canonical ensemble is 
the compressibility dnb/dfib and this is proportional to 
l/ab — {1 + a/N)/{2as) which is found to be positive 
even when N is set to unity. 

There have been several diagrammatic and fleld- 
theoretical approaches to the crossover problem [l^, [2^, 



42, 43, 44 



45|. Our results in Sections V are essentially 
the same as the diagrammatic approach of Hu, Liu and 
Drummond [Toj . although the derivations are somewhat 
different. In particular, in the diagrammatic approach 
the form of the gap equation was unchanged for con- 
venience and only the thermodynamic potential was al- 
tered. In our functional integral framework we can justify 
this as a natural approximation, and, in Section IX we go 
beyond this approximation and discuss the problems of 
self-consistently including the feedback of Gaussian fluc- 
tuations in the gap equation. 

The problems that we uncover in the self-consistent 
approach give insight into the conserving approximation 
scheme used by Haussman and collaborators They 
too flnd that, within their approach, as soon as one 
changes the gap equation from its mean-fleld-like form, 
one has problems with Goldstone's theorem. They fix 
this problem by simply redefining the scattering length in 
an ad-hoc manner to impose Goldstone's theorem. Our 
approach is fundamentally different, as it is based on the 
observation that the Goldstone mode is associated with 
the presence of soft phase modes in the superfluid phase, 
amenable to an amplitude-phase decomposition. 



XI. CONCLUSIONS 

To conclude, we have studied in this paper the BCS- 
BEC crossover in an attractive Fermi gas at T = which 
is relevant to experiments on ultracold gases with a wide 
Feshbach resonance. We have gone beyond the mean field 
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approximation and included the effects of quantum fluc- 
tuations at the Gaussian level. There is no small param- 
eter which controls this calculation, as we have not intro- 
duced a parameter like dimensionality (4 — e) or number 
of Fermion species 2N. Instead we have attempted to 
see whether there is an approximation scheme which can 
capture the known physics in both the BCS and BEG 
limits and in addition interpolate between them through 
unitarity. 

In summary: 

(1) Wc include the effect of quantum fluctuations which 
go beyond mean field theory using a functional integral 
approach at T = 0. We find that at the Gaussian level 
these fluctuations are the zero-point motion of the collec- 
tive modes and the virtual scattering of fcrmionic quasi- 
particles. 

(2) In the BGS limit, the virtual scattering of quasipar- 
ticles dominates the Gaussian correction and leads to 
Fermi-liquid corrections to the ground state properties. 

(3) In contrast, in the BEG limit the zero-point oscil- 
lations dominate the correction term. Wc can get an 
approximate understanding of the renormalization of the 
effective repulsion between molecular bosons and recover 
the Lee- Yang form for the quantum depletion. 

(4) At unitarity wc find that both collective modes and 
quasiparticle scattering contribute to the thermodynamic 
potential. Our results are in good agreement with both 
quantum Monte Garlo and experimental results. 

(5) We discuss in Section IX the problems of self- 
consistently including the feedback of fluctuations into 
the gap equation. Although the problem of impos- 
ing a gaplcss Goldstonc mode is solved by going to 
the amplitude-phase representation for the fluctuations, 
there are still some unsatisfactory aspects to the calcu- 
lation. One is the somewhat ad-hoc manner in which 
the ultraviolet divergences have to be regulated by inter- 
polating between the polar and cartesian representations. 
The results are not quantitatively superior to the simpler 
approach at unitarity and there is the further problem of 
thermodynamic instability in the BEG limit. In conclu- 
sion, we feel it is best to not modify the gap equation 
by feeding back the Gaussian fluctuations and to stick 
to the simpler set of equations dealt with in Sections V 
through VIII. 
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APPENDIX A: MEAN FIELD THEORY 

We review the derivation of the mean field gap and 
number equations with special attention to convergence 
factors, which will play a central role in a more compli- 
cated setting later on (see Appendix B). 

The saddle-point equation SSq/ 5 Aq = 0, with 5*0 given 
by cq. mi), leads to the MF gap equation 
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\i.,ikn 



dAo 



The Nambu Green's function 



Go(fc) 



1 



(^fc„)2 -El\ -Ao ^fc„ - Cfe 



(62) 



(63) 



is the inverse of Go ^(fc) defined in eq. p3p . Doing the 
Matsubara sum we get 

i.-iy I = (64) 
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To obtain the final result (|T4)) , we set the Fermi function 
/(-E-k) = at r = in the equation above, and use cq. ^ 
to take the infinite A limit and eliminate the coupling g 
in favor of the s-wave scattering length as- 

Evaluating n = —dflo/d^ with Qq = Sq/P we get 



(65) 



The Matsubara sum is formally divergent and we 
must introduce convergence factors. These factors arise 
because we need to calculate the equal time limit 
of G;i(k,r) = ^(rckT(r)4^(0)) and G02(k,r) = 

-(T4^(T)cki(0)) to obtain rikt = G?i(k,r 0') and 
rikj, = — G22(k, T 0+). We thus rewrite ((65)) as 



(66) 



as a contour integral in the complex z- 
z. The Fermi factors /(z) — XjieP^ + 

20 + 



and evaluate ^^j, 
plane, with ikn ^ 

1) ensure convergence for z — > +oo. For z — cx), e 
leads to the convergence of the first term but the second 
term is divergent. To convert the offending e"**"'"*^ in 
([Ml to the desired e'*^"" , we exploit the fact that the 
sum is over both positive and negative k and even under 
k -k, since G§2(-^') = -Gli{k) from eq. dM]). Thus 



1^ tanh(/3Sk/2) 

-c^k 



n = 2^G;,(fc)e^'="0^=^ 1 

]^,ikn k 

(67) 

The final result going from ([65]) to ([67]) could have been 
simply obtained by physical reasoning. The only point 
of going through the convergence factors in detail here is 
that it will streamline the discussion in Appendix B. 
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APPENDIX B: CONVERGENCE FACTORS FOR 
BOSE MATSUBARA SUMS 

In this Appendix, wc collect useful results for the 
asymptotic expansion of M^j for large frequencies 
and show that the integral of the phase S{q,uj) of 
DetM(q, + iO+) diverges for large negative uj. This 
forces us to introduce convergence factors to get a finite 
answer for the thermodynamic potential, leading us from 
the formal expression ([25]) to the convergent result ([26|) . 

We use eqs. (|21l22p to find the large qi expansion valid 
for qi ^ max(Ao, By neglecting the dependence on 
Ao and /i we get 



Mii{q,iqi) = A/22(q, -iqi) 
m \ - I 1 



47ra5 
m 



+ {m^^yAT:)^e^/2-iqi, 



2ek 



(68) 



where eq = q^/2m. For the slightly more restrictive case 
when qi further satisfies qi ^ eq/2, we get 



Mi2 = Alm^/\'/'/i2V27T). 



(69) 



On the real frequency axis, we are interested in large neg- 
ative w (at r = 0, the positive frequency contributions 
go to zero due to the Bose occupation factors nB{i^))- 
With |a;| » max(Ao, \fi\) wc find 



l3/2 

Anas 47r 
1- i- 

47r 



M22(q, w) 



Anas 



771 ' / 



(70) 



Here the imaginary part of Mn comes from just the first 
term in the sum in (|2ip . for which the branch cut appears 
on the negative real frequency line. Using the slightly 
more restrictive condition |cij| ^ eq/2 

Ml2(q,c.) = -^2!!^|c.|-3/2(l + ^) (71) 
in 

Thus in the limit of large and negative uj we can neglect 
Mi2 and get to the leading order 



DetM 



777^/^ r 777,'^ 



16n^as 



16n^ 



(72) 



The phase is given by Im In DetM = (5 « 
tan~"'"(a5-\/77i|(jj|) « n/2 as w ^ — oo. Thus the Mat- 
subara sum in eq. ()25p is divergent. 

We next turn to the derivation of the convergence fac- 
tors in eq. (pS)) and how they lead to finite results. We 
begin with looking at the sum on Matsubara frequencies 



for a fixed q: 

^ In DetM (q) 



E 



In Mil + lnM22 + In 1 - 



M?2 



M11M22 



(73) 



and show that the first two terms should be written 
as J2iqi flnMiie**°^ + In M22e-"""^ 1 , while the third 
term does not need a convergence factor. We can rewrite 
eq. (IHI) as 

Mil (q) = M22 i-q) = - + E (74) 

with GO = GO(fc) and G° ' = G°{k + q). We then ex- 
pand the logarithm in powers of g so that In M22 — 
Eg [ln(l/g) + 9 Efc G§2G0'i + ...]. Using the argument 
given in Appendix A below eq. ([55)) . we see that the equal 
time limit requires that G22 carries a factor of e"'"" 
and G'li a factor of 6'^'^"+'''^° . We thus see that or- 
der by order in g each term in In M22 comes with a fac- 
tor of e~*'"° and InMn comes with e+**'° . We also 
note that, using this prescription, Mi2(<z) = 'M.2i{q) = 
Ek ifc G52G52 does not acquire a convergence factor 
and, in fact, none is needed. 

The Matsubara sum Eiq, is converted to a standard 
contour integral. Convergence for z — s- -|-oo is guaranteed 
by the Bose function nB{z) = l/{e^^ — 1). For z —00, 
convergence is ensured by converting the problematical 
factor of e"*'''" into the convergence factor 6*''° , fol- 
lowing the same reasoning as in Appendix A. Using the 
fact that the sum is over both positive and negative q 
and M22(<z) = M.ii{—q) (see eq. (|2T|) ) we obtain 



^ In DetM (g) 

iqi 

= In 



21nMiie^«'°^+ln(l- 



M22M11 



MiiDetM' 



M. 



(75) 



which is exactly the expression in eq. ()26p . 

We finally show explicitly that the Matsubara 
sum in eq. (|75|) is convergent. The Matsubara sum 
in eq. (|75|) can be written as the contour integral 

dz/{2ni)nBiz) ln[Mii(q, 2;)DetM(q, z)/M22(q, z)] 
where C runs on either side of the imaginary z axis, 
enclosing it counterclockwise. We distort the contour to 
run above and below the real axis and at T = obtain 
for the thermodynamic potential 

^9="^^E / — ['^(q,w) + <5ii(q,w) - (522(q,t^)] 

(76) 
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where (5(q, t^;) = ImlnDetM(q, w + iO"'") and Sn and S22 
are the corresponding phases for Mn and M22- 

From the leading order expression for Mn in eq. ([70|) . 
we see that for large negative lo^ Su 



and hence 

that term is convergent. To look at (5 — 622, we recognize 
that this is the phase of DetM/M22 = Mn - Mf2/M22. 



1-7/2 



Now for large negative iv, 'M.I2/M22 ^ 
and can be neglected in comparison to Mn. Thus the 
integrand in eq. ((76|) reduces to 26ii and we get a con- 
vergent answer. 



APPENDIX C: NUMERICAL EVALUATION OF 
BOSE MATSUBARA SUMS 



While the real frequency representation of cq. (j26|) 
gives physical insight into the deviations away from mean 
field theory, it is numerically simpler to do the calcula- 
tion on the imaginary frequency axis. On the real axis 
one encounters principal part singularities analogous to 
the ones encountered, e.g., in the normal state calcula- 
tions of ref. [4^ but further complicated by the broken 
symmetry in the supcrfluid state. 

If one wants to use eq. on the Matsubara axis, 
one needs to explicitly take into account the convergence 
factor e"'"*''''^ and take the limit r — > 0"*" at the end. Here 
we outline an alternative procedure which simplifies the 
numerics. Let us begin by looking at a part of Mii(q, z) 



Mf,(<z) = i + E 



E-E' 



which has no singularities (poles, branch cuts) or zeros 
in the left-half plane (Rez < 0). Since we will use this 
to get convergent results wc call it Mfj^ and M^2- 
We may write the Mn piece of eq. (|75p as 



E21nMiie 



E 



2[ln(Mii/Mfi)-f InMfi] 



(77) 

where we drop e'*'° on the right because each term is 
convergent. The Matsubara sum of the second term is 
seen to be zero at T = by evaluating it as a stan- 
dard contour integral and noting that In Mfi has no sin- 
gularities in the left-half plane. In fact, now we may 
write the above result in a more symmetrical form as 
X;^,, [ln(M22/M^2) + ln(Mii/Mfi)] and combine this 
with the second term of ((75|) to obtain 



E In DetM((7) E 1^ Det 



M{q) 



MC{q) 



(78) 



where the matrix M*^(g) is a diagonal matrix with the 
entries M^2('Z) ^^d 'M.f^{q). This expression leads to a 
rapidly convergent answer, which in the T = limit can 
be evaluated as an integral along the imaginary axis in 
the iqi ^ z = (x + iy) plane with f3~^ Sig, ~^ I dy / {'^tt) . 



+ 



\/+- 



Hartree 
term 



n th term 



FIG. 10: Diagrammatic representation of the Gaussian cor- 
rections to thermodynamic potential in the BCS limit. The 
full lines are fermion propagators and the wave lines represent 
the attractive interaction. The first diagram corresponds to 
the Hartree term. 



APPENDIX D: FERMI LIQUID CORRECTIONS 

FOR kplas] < 1 WITH as < 



Here we give some details of the argument that shows 
that the Gaussian corrections to the thermodynamic po- 
tential in the extreme BCS limit of the attractive Fermi 
gas have the same expression as_the _standard Galitskii 

[ii. 



25, 2 



of the repulsive 
First we recall 



and Huang-Lee- Yang theory 
Fermi gases with a sign change in as 
the usual Galitskii theory and discuss why it is not di- 
rectly useful for as < 0. Next, we describe how the BCS 
limit results of the superfluid state theory developed in 
the paper are related to those of normal state Galitskii 
theory. 

In the Galitskii theory of the normal Fermi gas the 
thermodynamic potential f2 is written in terms of the 
two-particle propagator, which is the sum of particle- 
particle channel ladder-diagrams of Fig. [101 For a re- 
pulsive interaction Vq, wc find 



f^.ce + 2E(-l/^E 

1=1 q 



J2Goi~k)Go{k + q) 



(79) 

where Gq is the non-interacting Green's function Go(fc) = 
(ifc„ — Ek + m)~^- To make contact with results of our 
paper, it is useful to sum up the series and write it as 



^hee + 2 E - E Go{-k}Go{k + q)]. (80) 



The repulsive Vq can then be replaced by as > in the 
usual way using m/AiTas ~ l/Vg -I- l/^^q- 

It is well-known that the pairing instability of the nor- 
mal Fermi gas to attractive interactions implies that we 
cannot extend the Galitskii calculation directly to the 
case of attractive interactions < 0. If we were to try 
and set Vq = — 5, the attraction of eq. ([2]), we would 
find that, for small g, there is a pole on the imagi- 
nary axis in the upper-half plane in addition to a branch 
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q)] , and wc have used convergence factor tricks to obtain 
well-defined Matsubara sums. The last expression is the 
same as the phase shift obtained from the normal state 
result above ([50)1 with Vq ^ —g. We thus conclude that 
in the Aq — > hmit of the superfluid state, the Gaussian 
correction to H, is the same as the Galitskii and Huang- 
Lee- Yang result for the repulsive Fermi gas with sign of 
as changed to as < 0. 



APPENDIX E: THERMODYNAMIC POTENTIAL 
IN THE BEC LIMIT 



FIG. 11: Analytic structure of two-particle propagators for 
the attractive Fermi gas. The upper panel correspond to the 
unstable normal state which has a branch cut for to > —2^ -|- 
eq/2, representing the continuum of excitations, and a pole 
on the imaginary axis in the upper-half plane, for small q, 
which signals the BCS pairing instability. The lower panel 
shows the analytic structure of M~^ in the stable superfluid 
state with poles at ±t<;o(q), the collective mode frequencies, 
and branch cuts beginning at ±min(_Ek -I- £^k+q). 



cut along the real axis ui > — 2yU -f- eq/2, as shown in 
the top panel of Fig. [TlJ This pole, which occurs at 
z ^ +iefexp{—l/kp\as\) for q = 0, is the signature of 
the Cooper pairing instability. 

The superfluid state calculation for the attractive 
Fermi gas discussed in the text of this paper deals with 
the broken symmetry saddle point which is stable. The 
two particle propagator in the superfluid state has a 
very different analytical structure (lower panel of Fig. [TT|) 
compared with the unstable normal state just discussed 
(upper panel of Fig. [11]). As described in Section IV be- 
low eq. (PH]) . M~-^ has poles at the collective modes fre- 
quencies and branch cuts corresponding to the gapped 
two-particle continuum. 

In the BCS limit, Aq ^ and the contribution of the 
zero-point motion of the collective modes to f2 is utterly 
negligible for phase-space reasons, as discussed in Sec- 
tion VI. Thus the two-particle continuum dominates the 
Gaussian correction tofl. In the limit Aq the branch 
cut extends over the entire real line, although the appro- 
priate limit of the itfc, Vk factors shows that the phase shift 
vanishes for uj < —2/^ + eq/2. In this sense the contin- 
uum contribution to ft given by J_^''^'^'' dujS{q,uj) may 
be simplified with the lower limit becoming — 2/i -f- eq/2 
and the upper limit going to zero. 

Now one can check that this continuum contribution 
is exactly the same as the corresponding continuum con- 
tribution of the normal state Galitskii theory, changing 
the sign of as- Note that the singular pole piece does not 
show up in this result. To see that the (5(q, w)'s in the 
two theories are the same, start with the M matrix of the 
broken symmetry theory. As Aq — > 0, M12 vanishes and 
<5(q,w) ~ 2ImlnMii = 2Imln[l - g^^fc Go(-fc)Go(fc 



Here we bricfiy sketch how we arrive at the leading 
order correction to the thermodynamic potential in the 
BEC limit eq. ([55]) starting with the results of Appendix 
C. In the BEC limit, fi < and approaches half the 
binding energy of the molecules = l/(2ma|). Thus 
Ao/|/i| <C 1 and can be used as an expansion parameter. 
Then one can write = 1 — Aq/4^^, — Aq/4^^ and 
-£'k = + AQ/2^k with = ek l/^l to leading order in 

Ao/ImI- 

Now writing Mn = Mf^^+SMn, where Mf^ is defined 
in eq. ([77)) and 



<5Mn(g) 



E- 



iqi+E + E' 



(81) 



one can easily see that ^Mn ^ Ag while M12 ~ Ag. 
Then, to order Aq, the expression in cq. ([78]) gives 



^Re 



M?2 



SMu 



(82) 



To leading order, SMu = A^/\n\y^F{q/ y/2m\fi\,qi/\fi\) 
where -F(Q, Qi) is a dimensionless function given by 

^ (83) 



16 ^ iQi 



2 + 2A-2 + Q2/2 
1 



((K + Q/2)2 + l)2((K-Q/2)2 + l)2 

Here we use capital letters for dimensionless vari- 
ables Q = q/^y2m\f^\ and Qi = q;/|/i|. M12 = 

A'a/\f^\^^^Ii'i/ V^MiA^ «/ImI) where /(Q, Qi) is a dimen- 
sionless function given by 

KOO)- 2 + 2i^^ + QV2 

1 



((K + Q/2)2 + l)((K-Q/2)2 + l) 
and Mfi = vM^(q/VM,9i/lMl) where H{Q,Qi) is a 



dimensionless function given by 



1 



H{Q,Qi) = — V-*Qi/2 + OV4 + l--= 7= 

on \ \/2mas\/\l-L\ 



^ (V-*Qi/2 + QV4 + i-i) 



(85) 
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Putting all these together, we get 



{2m) 2 

47r3|/i|3/2 



d^Q / dQ 



Re 



F{Q,Qi 



H{Q,Qi) 



2\HiQ,Qi)\' 



(86) 



Numerical evaluation of the integral gives f2g = 
-a(2TO)3/2A;^/|^i|3/2/(2567r) with a = 2.61. 



APPENDIX F: THE AMPLITUDE-PHASE 
ACTION 

Starting with A(a;) Ao[l + A(a;)]e**(^) (see eq. (P])) 
we transform to a gauge where A(x) is real. We trans- 
form the fermion fields 



(87) 



with 



U(x) 



-ie(a;)/2 Q 





so that the action ([7]) now reads 



1 



•^^.A = / d2;-|A(a:)|2-^t(3,)G-i(a:,x')V(a;') (89) 



where G"i = UGUt. ^ 

We can now write G^^ = Gq ^ + K, where G,^^ is 
the (inverse) Nambu Green's function defined by ([8|) in 
k ~ (k, ifc„)-space. The matrix K is 



K{x,x') = 



AoX{x)a, + ^ (w0{x) ■ V + ^W^d{x) 



Ldre{x) + ^\V0ix)\']a, 
2 8m 



6{x ~ x') (90) 



whose Fourier transform is 



K{k',k) 



S{k -k' + q) 



+ 8to ^ ' 92^*91 '^3 (5 (fc - fc' - gi - ga) (91) 



9192 



with (5^ = ^k+q - Ck- 

Integrating out the fermion fields -0 we get the func- 
tional integral (j50p with the action SAo,x,d = + Sg of 
ea ([5T|) . The 5*0 piece, defined in (|12p . comes from the 
Gq ^ term; for the J term in (|50p see ([5^. The Gaussian 
piece, arising from K, is given by 

Sg = ^Xg\-g - TrGo(A:)K(fc, k) 

+ irrGo(fc)K(fc, A: + g)Go(fc + g)K(fc + q, fc) (92) 



The Gaussian action of (|53p follows immediately from 
(|9T|) and (l92|), with the D matrix given by dH]). Our 
next task is to derive the equivalent expression for the 
D matrix (|55p which is written purely in terms of G*^, 
without any iqia^ — 6^ factors. The case of Dn is sim- 
ple; there are no such factors to begin with and we only 
need to evaluate the Nambu trace in ([54]) to obtain ([55]) . 
In what follows, we use the notation G = Go(fc) and 
G' = Go(fc + g), and drop the subscript for notational 
convenience. 

In order to write D12 and D22 in terms of the Green's 
functions, one needs to express the vertex iqia^, — (5^ in 
terms of matrix elements of G~^. It is easy to see that 
the vertex can be written as 

■^11 

-(522 



iqicrs - (5^ = V = 



G22 . We will 



where Su = G[^^ — G-^^^ and 822 = G'22^ 
also use the following identities: 

G22G22^ = 1 — A0G12 , GiiGj^j^ = 1 — A0G12 

Gi2G22^ — — AqGii , Gi2Gj^j^ = — A0G22 (93) 
By definition we have 

D12 = Ao^TrGo(fc + g)VGo(fc)ai 

= ^ 5ii(GiiG'j2 + Gi2G'^j^) 

fc 

-<522(Gi2G^2+G22G'i2) (94) 

Now using the identities of eq. ((93)) . we get 
J2k '^ii(GiiGi2-fGi2G'ii) = J2k Ao(G22Gii — G11G22) 
and a similar result holds for the 822 piece. Adding both 
terms we get D12 of eq. (|55|) . For D22, one can write 



z 1 
D22 = |- V(Gn - G22) + -TrGo(fc)VGo(fc + q)V 

k 

(95) 

The second term above can be written as 

\ Ylk <^iiGiiG'ij+(5|2G22G22— 2(5ii(5iiGi2Gi2 
Using the identities in cq. ([M)) . we get 

GiiG']^]^(52i = A5 G11G22 + G']^jG22 — 2Gi2G']^2 

k k 

+2Ao^(Gi2 + G;2) 

k 

-l-^GiiG'jj^ + G'liGjj^^ - 2 

and a similar result holds for the 622 piece. For the last 
term we get 

— 2 Gi2G']^2'^11'522 = 2Aq [G11G22 + Gj]^G22 

fc k 

-2G12GI2] + 2Ao ^(Gi2 + G'la) 
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Adding all the terms we get 

^ G22G'ii + G22G11 — 2Gi2G'^2 "I — 2" $Z 

k k 
+ - ^ G22G22^ + G22G22"^ + GuG'^j^"' + G^^Gj/ - 4 



where we have used G12 = X^fc ^'12 • We now use 
G'l^i = Gn' + ig, - 5£, and G'22' = G22' + iqi + 
to write the last term as —■^^oj^k^i^ + ■^"(Gii + 
G22) - ^{Gii - G22 - G'li + G22), where we have 
used J2k = Sfc G'ij to write this form. Now us- 
ing proper convergence factors Gn + G22 = and 
so the terms multiplying iqi vanishes . In the last term, 
replace k ^ k + q to show that this term is propor- 
tional to q^/2m and it actually exactly cancels the sim- 
ilar term in D22 coming from the (q^/2TO)TrGcr3 piece. 
Now J2k *^i2 = — Aq J2k DetG, and so combining every- 
thing we get the result for D22 in ([55)) . 

D22 ~ y ] G11G22 + G']^]^G22 — 2Gi2G^2 ^ 2DetG 

(96) 

Going to the rescaled basis (A, 9), we then have 

Dii = — h - [G22G']^]^ + G11G22 -I- 2Gi2G']^2] 
■9 ^k 

622 = ^^ [G22G'ii + G^2Gii - 2Gi2G'i2 - 2DetG] 



D12 — [G22G'jj - G11G2; 



(97) 



Just as in the case of the static Saddle Point number 
equation, one runs into formally divergent quantities in 
evaluating the q sum to get the action. To fix these, one 
has to regularize using proper convergence factors. The 
D basis is not the basis of choice for fixing the conver- 
gence factors. Instead of the amplitude A and the phase 
9 we can work with the complex fluctuation fields defined 
in (1561). In this basis the matrix D is transformed to 



L = WDWt 



(98) 



Since this is an unitary transform DetD = DetL. Then 
we get 



L12 — Gi2G']^2 



X 


= Mn- 


X 






^ y 




X 


= M12 + 


X 


+ y 








X 


= M22 - 


X 






^ y 


T 



(99) 
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where X = 1/g + J2k DetG = 1/g - J2k l/(2£^k) is the 
LHS of the mean field gap equation and L12 = L21 . 



Now we can fix convergence factors with InLn carry- 
ing a convergence factor of 6"*"'''° and lnL22 carrying a 
convergence factor of 6^*"^'° . This can be seen by ex- 
panding the log and remembering G22 carries a factor 
of 6^'*^"° , G[i carries a factor of 6+'^''"+'")'' and so 
on. The reasons for the convergence factors are related 
to taking the correct equal time limit and is discussed in 
in detail in Appendices A and B. 

Remembering Lii(— g) = L22('7) one can take out 
In Lii +ln L22 from the In DetL and then convert In L22 to 
InLii using q —q. Now one can convert the Matsub- 
ara sums to real frequency integrals which are convergent. 
The resulting action is 



Sn 



Lii 



DetL 



(100) 



One can then follow the asymptotic forms of the M ma- 
trix derived in Appendix E to get the large energy-short 
wavelength behavior of the L matrix. The asymptotic 
forms of Lii and L22 are the same as that of Mn and 
M22, with m/Anas replaced by m/A-Kas — X/2. How- 
ever the presence of the X/2 factor in L12 fundamentally 
changes its asymptotic behavior from o;"'^/^ for the M 
matrix to a constant {X/2). We can thus no longer ne- 
glect the L12 terms in the high frequency limit and this 
leads to a divergent answer. 

APPENDIX G: DILUTE BOSE GAS 

In this appendix we show how the method of Gaussian 
fluctuations yields the correct answers in a somewhat dif- 
ferent problem, that of a dilute Base gas with repulsive 
interactions. Although the results are standard [26|, the 
method used here parallels that used in our paper, and 
serves to illustrate several technical points including: (1) 
the role of convergence factors, (2) retaining the mean 
field form of the saddle point equation and including 
quantum fluctuations in the thermodynamic potential, 
and (3) taking into account the ^-dependence of the the 
saddle point in the number density equation. 

The Hamiltonian for a repulsive (go > 0) Bose gas is 

H = J d^x<i>*{x){-V^/2M ^ fx)'^{x) 



-hy$*(a;)$*(x)4'(a;)$(x) 



(101) 



Writing <I>(.t) = $0 + (i^) the action is S ^ So + Sg- 
where 

So = /? + 
and the Gaussian part is given by 

^. = |E(C*('z),C(-'z))A(g)^ ^('^^ 



Ci-q) 



(102) 



(103) 
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Here Aii(g) = A22{-q) = 
eq = |q|V2M, and Ai2(<z) 
ing out the C fields we get the thermodynamic potential 



-m + Cq - 

A2l(-(?) 



H + 250*0' with 
= go*o- Integrat- 



~ f^o + (1/2/3) InDetA(q), 



(104) 



where f^o = So/p. The Matsubara sum in the Gaus- 
sian piece is ill-defined. We write lnDetA((7) = In An + 
lnA22 + In (l — Af2/-^22Aii), introduce convergence 
factors of exp(ig/0"'') with the An term (associated with 
C*C)j and exp(— i(7;0+) with the A22 term (corresponding 
to CC*) and use q ^ —q to write the A22 piece in terms 
of All. At T = the sum Xl^,, 2 In An exp(ig,0+) 
vanishes by contour integral methods since the integrand 
has no singularities in the left-half plane. The remaining 
sum can be explicitly done by contour methods to obtain 



2.9o$g) 



(105) 



where i?q = ^(eq - n + 2go^l)'^ - g^^"^ is the Bogoli- 
ubov dispersion. The quantum fluctuations arc clearly 
seen to have the form of zero point motion of the collec- 
tive modes Eq/2 with a "convergence factor" subtraction 
which eliminates the ultraviolet divergence by canceling 
out the contribution of the quadratic part of the Bogoli- 
ubov spectrum at large q. 

The uniform, static saddle point is determined by 
6So/5^o = 0, so that 



We now solve this equation for /Lt(n) in powers of {nal). 



To leading order fi 
Of, we get 



Airnab/M and to the next order in 



fj, = 47rnab/A/[l 



-32/37r-i/2(ria3)i/2] 



(109) 



This is the correct equation of state for a Bogoliubov 
dilute Bose gas, including the Lee- Yang correction. 

We note that we cannot identify the saddle point value 
of $0 in eq. (|106[) with the condensate fraction, once 
quantum fluctuations are taken into account. This iden- 
tification is usually made, together with the replacement 
Po ~> 47raf,/M. We note here that this identification 
makes (|105p divergent and is thus not well defined.. To 
find the condensate fraction we use the expression for the 
q 7^ momentum distribution 



n(q) 



to derive the well known result 



(A-^)ii 



(q^O). 



(110) 



(111) 



From the quantum depletion we can obtain the well 
known result for the condensate fraction using iVg = 

^ - Eq^o "(q) = ^ (1 - 8/3v/(na,)V^) . 



and this condition is again needed in order to satisfy that 
the excitation spectrum is gapless. We use {d^l/dy) = 
—N to determine /x. In evaluating the thermodynamic 
derivative we cannot treat $0 as a constant, and must 
keep track of the //-dependence of $0 hi eq. (jl06p . We 
thus get n = -[i^ jlgQ + 5 Eq (-^'q - ^q ^ m) with = 
\/ (cq + /^)^ ~ M^- Taking the derivative with respect to 
/i, we get n = fj,/go + \ Y^qi^ ^ Cq/^q)- 

Now, using the relation between the bare repulsion g 
and the boson scattering length at, given by M/Airai, = 
1/50 + Eq l/2eq we get 



i^/Eq - 1/eq) 



where we have added and subtracted /i/2_Eq to isolate 
the cancellation of divergences. The first integral is 
-(l/37r2)M3/2;x3/2 and the second one is -Af3/2^3/2/^2^ 

So in all we get 



n — fiM/4:T:af, 



(4/3^2)^,^3/2^3/2 
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